Revealing the spatiotemporal complexity of the magnitude distribution and b-value during an earthquake sequence

The Magnitude–Frequency-Distribution (MFD) of earthquakes is typically modeled with the (tapered) Gutenberg–Richter relation. The main parameter of this relation, the b-value, controls the relative rate of small and large earthquakes. Resolving spatiotemporal variations of the b-value is critical to understanding the earthquake occurrence process and improving earthquake forecasting. However, this variation is not well understood. Here we present remarkable MFD variability during the complex 2016/17 central Italy sequence using a high-resolution earthquake catalog. Isolating seismically active volumes (‘clusters’) reveals that the MFD differed in nearby clusters, varied or remained constant in time depending on the cluster, and increased in b-value in the cluster where the largest earthquake eventually occurred. These findings suggest that the fault system’s heterogeneity and complexity influence the MFD. Our findings raise the question “b-value of what?”: interpreting and using MFD variability needs a spatiotemporal scale that is physically meaningful, like the one proposed here.


Supplementary Notes Note 1: Cluster size statistics
Supplementary Fig. 2 summarizes the cluster statistics in terms of size and relative proportions for each period, not accounting for short-term aftershock incompleteness (STAI, Note 2.4).
From absolute cluster sizes ( Supplementary Fig. 2a), it is evident that • the 'Norcia' is the most active period; • only two or three clusters are simultaneously active during the 'Amatrice', 'Visso', and 'Norcia' periods • all clusters are about equally active in the 'Campotosto' period; The temporal distribution of each cluster ( Supplementary Fig. 2b) reveals that: • after the Norcia event, all clusters experience a "boost" in seismicity; most of the earthquakes inside and outside the clusters belong to this period, except for C4 and C5; • C4 and C5 are most active in the Campotosto period.
The relative cluster size ( Supplementary Fig. 2c) shows that: • in each period during the sequence, ≳ 50% of the earthquakes belong to a cluster; • C1 is the only cluster that is constantly active and never quiet.

Note 2: Additional aspects of MFD and -value analysis
With particular reference to our analysis, we here discuss further aspects that we consider important when performing and interpreting MFD and -value analyses.

Sample size and magnitude range
First and foremost, robust -value estimates with low uncertainty require a sufficient sample size. In fact, the maximum likelihood estimation assumes that the GR relation holds over a wide magnitude range, typically at least three magnitude units 2 . A MFD that spans less than this range biases the -value estimate due to a correlation with the largest magnitude and, related to it, a truncated or tapered MFD 3 . In our analysis, we are close to this limit but could not reduce the overall cutoff magnitude below w 1.5 because the Lilliefors c over the whole sequence is not below this level, even in the complete periods, i. e., those unaffected by STAI (see Note 2.4). Although we did not model the MFD with a tapered GR distribution, we can exclude the influence of the corner magnitude on -value variations because we do not observe a dramatically varying maximum magnitude among the clusters or their temporal subsets. MFD and -value variations therefore originate from differences or changes in the entire exponential part of the MFD.

Exponentiality, Lilliefors test, and -value stability
For a -value of 1.0, a magnitude range of three units corresponds to a sample size of 1000. For such rather large sample sizes, the Lilliefors test may still underestimate the departure from exponentiality, i. e., not detect those even if a departure exists 3 . Inferring the robustness of the -value estimate therefore requires a careful analysis of both the MFD exponentiality and the -value behavior as function of c .
Lilliefors c may only relate to a short-lived exponentiality (i. e., no persistent exponentiality with increasing c ) 4 , for instance in the 'pre-Visso' and 'pre-Campotosto' period of C3 (see Figs. 4 and 5): at a low magnitude level, the -value indicates an exponential behavior ( ≥ 0.1), but a departure of exponentiality ( < 0.1) at a higher magnitude level. Such short-lived exponentiality may either be due to random fluctuations (especially for small sample sizes) or may have real physical or technical causes (e. g., scaling breaks, incompleteness, or magnitude artifacts) 3,4 . In other cases, Lilliefors c indicates an exponential behavior, but the -value continuously increases with c (e. g., for the 'mid' period, see Supplementary Fig. 7)-this behavior is an indication for STAI (see Note 2.4). For this reason, an inspection of the -value stability with c is always advised.

Influence of magnitude scale
For a consistent analysis of the MFD, moment magnitudes, w , are preferred over local magnitudes, L , as discussed in Herrmann and Marzocchi [4].
Supplementary Fig. 3 shows that using L introduces a different MFD behavior for the individual clusters. Compared to w -based MFDs (Fig. 3), exponentiality of C1 and C2 is rejected for a wider magnitude range. For C3, the exponentiality is only short-lived below L 2.0. This behavior is explained with a stronger under-representation of the L scale 4 , likely due its intrinsic scaling break toward low magnitudes 5,6 . The -value at Lilliefors c differs not only in a relative sense (i.e., among the clusters), but also in an absolute sense (being much lower).
Using L for the individual periods, this larger discrepancy is reduced (see Supplementary Figs. 4 and 5) and reproduces our main findings that we obtained for w (i. e., that the MFD can vary spatially and temporally, but can also remain constant, and that the b-value increases in C1 prior to the Norcia event). However, the -value behaves differently as function of c , owing to the scale change: For C1, the -value does not increase as strongly anymore with decreasing c ; for C2 and C3, the -value is not constant anymore with c , but decreases with decreasing c . w are important to guarantee consistent analyses such as to obtain a more robust indication for departure of exponentiality (by excluding departures due to L 's scaling break). Yet, it remains unclear if the magnitude conversion relation of Grünthal et al. [7] used by Tan et al. [1] for converting L into w represents a suitable choice; it will eventually saturate toward low magnitudes (i. e., overestimate w ) owing to its polynomial form. This saturation may explain the apparent over-representation at low w , favoring an increasing -value with decreasing c as observed in many subsets and in the overall catalog with STAI periods removed ( Supplementary Fig. 9).
The relation of Grünthal et al. [7] may not be ideal, but we wanted to use the simplest solution: taking the original w estimates provided with the catalog. Alternatives are the relation used in Lolli et al. [8], i. e., the HORUS catalog, a revision of Gasperini et al. [9], or the bilinear relation of Malagnini and Munafò [10]; but using either of those would increase the entropy of our analysis and inevitably the "influence of expert choices" as stated in the main text-a condition that we want to avoid as much as possible as advocated by our study. In this sense, the original catalog creators did the choice for us. The relation of Grünthal et al. [7] is not unsuitable: it is based on European data, which includes data of Italy . Tan et al. [1] did not use this relation directly, but calibrated the constant (using ∼500 earthquakes with w estimated from regional waveform fitting), as specified in their supplemental material.

Short-term aftershock incompleteness (STAI)
It is well known that STAI has a large influence on the perceived -value change 11,12,13 . The Norcia mainshock has the strongest effect, see Supplementary Fig. 7 ('mid' period): For all clusters, the -value in the 'mid' period continuously increases with increasing c , indicating a gradually curved MFD due to incompleteness. Recall from Fig. 4 that the periods excluding STAI had shown the opposite behavior, highlighting the strong influence of STAI on the MFD shape by reducing the estimated -value. In C1 and C2, the 'early' period has lower -values than their 'pre-Visso' period ( Supplementary Fig. 7), which purposely excluded the STAI effect of the Amatrice event; in C3, this effect is not evident, indicating that the Amatrice event is too far south to influence the detection of earthquakes in the northern part. The Campotosto events in the 'late' period also cause a STAI effect (see Supplementary Fig. 8), albeit only in C2, which is located closest to the Campotosto events in the southern part; C1 and C3 are too far away to observe the same effect. Given these observations, STAI needs to be carefully addressed, e. g., by disregarding incomplete periods, but only as much as necessary to not reduce the number of usable earthquakes unnecessarily. An alternative is using only positive magnitude differences for estimating the -value because they are not affected by STAI 13 .

Dependence on the used catalog
In Supplementary Fig. 9, we compared the w -based magnitude statistics of Tan et al. [1]'s high-resolution catalog with the HORUS catalog for the CI2016 sequence. Event selection comprises the spatiotemporal extent of the extracted high-resolution catalog, i. e., from 2016-08-15 until 2017-08-15 within the volume it covers. The comparison was performed for two cases: (i) using all events during the sequence, and (ii) using only events in non-STAI periods that are indicated in Supplementary Fig. 6. The analysis shows that the -value differs considerably between both catalogs (0.2 units at Lilliefors c ) despite using the same magnitude types ( w ). Therefore, the w estimates are not consistent with each other. Incompleteness may be ruled out as a cause, since the -value also differs by the same amount in the non-STAI periods. Due to this inconsistency, estimated -values cannot be compared among these catalogs, for instance to relate the estimated -value during the sequence using the high-resolution catalog with a reference "background" -value using HORUS. This finding highlights the need to treat w obtained with different conversion relations with utmost care in statistical analysis.
More generally, the choice of the used catalog may influence the analysis because catalog creators often use different approaches to determine earthquake locations and magnitudes.

Possible contamination by quarry blasts and explosions
An earthquake catalog may be contaminated by unrelated events. Based on recent findings and a discussion that quarry blasts contaminate the MFD and bias the -value 14,15,16 , we obtained locations of quarry blasts in the surrounding of the sequence since May 2012 until August 2021 (see Supplementary Fig. 11). No identified quarry blast or explosion occurred during and within the extracted catalog of the CI2016 sequence. Yet, for good measure, we repeated the analysis shown in Fig. 4 using only night-time seismicity (19:00 to 05:00 UTC, see Supplementary  Fig. 13), but cannot find an impact on the -value trends discussed in the main text. for hypocenter density (within a cubic kilometer, km -3 ; see color bar). The density is shown at every hypocenter whose corresponding density is higher than 20 km -3 , which is the median density among all hypocenters. The depth sections are to scale. The main events 'Amatrice', 'Visso', 'Norcia', and four times 'Campotosto' are represented by larger blue circles, annotated with the respective initial letter (A, V, N, C) in the map view.     [17]. The yellow-red overlay represents a kernel density estimate of the shown points, which reveals the decay of the short-term aftershock incompleteness (STAI) effect to a background level (suggested by W. Ellsworth, pers. comm., 2021). This visualization helps us to define the necessary "safety margin" to reduce the bias due to incompleteness by only using events in the gray-shaded periods (individually for each period).   ingv.it). The performed analysis is analog to Fig. 3 and Supplementary Fig. 3. For both catalogs, the analysis is performed for (i) all events during the sequence, and (ii) events during the sequence without the periods affected by STAI after each main event (i. e., only using the gray-shaded periods in Supplementary Fig. 6), including the main events themselves. See Supplementary Note 2.5 for more information and a discussion of these results.  Fig. 1, but showing only the temporal subsets before and after the Visso event ('pre-Visso' and 'pre-Norcia', respectively) for Cluster 1-3 (see legend). Note that the geographical bounds slightly changed compared to Fig. 1 by zooming in. The depth sections are to scale. The figure shows that (i) Norcia nucleated in between those two subclusters of C1; and (ii) C1's pre-Norcia seismicity occurred mostly at depth, along the north-eastern extension of the subhorizontal shear zone.